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Abstract 

We present a method to compute high-order derivatives of the total energy 
which can be used in the framework of density functional theory. We provide 
a proof of the 2n + 1 theorem for a general class of energy functionals in 
which the orbitals are not constrained to be orthonormal. Furthermore, by 
combining this result with a recently introduced Wannier-like representation 
of the electronic orbitals, we find expressions for the static linear and nonlin- 
ear susceptibilities which are much simpler than those obtained by standard 
perturbative expansions. We test numerically the validity of our approach 

with a ID model Hamiltonian. 
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Perturbative techniques are usually applied to density functional theory (DFT)Bi to 
study the response properties of materials from first principles. The evaluation of the second- 
order derivatives of the total energy yields phonon spectra,! effective charges,!'! dielectric 
constants,!'! piezoelectric tensors!! and many other experimentally measurable quantities. 
Likewise the computation of higher-order derivatives permits the ab-initio prediction of 
properties such as the Raman tensors,! the second and higher-order susceptibilities, the 
nonlinear elastic constants etc.. 

There are already very sophisticated analytical methods to obtain the values of the 
second-order derivatives, and today it is possible to evaluate these quantities in systems 
with many atoms per unit cell.0 On the contrary, the evaluation of third or higher-order 
derivatives relies mainly on finite differences: the required derivatives are computed by nu- 
merical differentiation of the second-order derivatives. The cost of the finite differentiation 
limits the applicability of the technique to small systems and to short wavelength perturba- 
tions. 

Closed-form expressions of the third or higher-order derivatives, obtained by a straight- 
forward application of quantum-mechanical perturbation theory, are usually cumbersome. 

In the case of the second-order susceptibilities, i.e third order derivatives of the energy 
with respect to a uniform electric field, the perturbative expansion provides a formula which 
apparently diverges in the static limit. This divergence can be eliminated as shown in 
Ref. ( [TTJ) in the context of a non self-consistent electronic structure theory. A specific 
application of the resulting formula has been performed using a semi-empirical tight-binding 
Hamiltonian.0 In order to extend this scheme to self-consistent DFT one has to face rather 
formidable formal difficulties. An explicit expression for the second-order susceptibility 
within DFT has been obtained using a software package for symbolic manipulation^ and 
it has been applied by just one research group due to its complexity!! 

Alternative analytical expressions for the high-order derivatives of the energy are pro- 
vided by the 2n + 1 theorem, well known in standard perturbation theory for many years0 
and recently rewritten in the language of DFT.0 This theorem states that the derivatives of 



the energy up to order 2n + 1 can be computed if the change of the wavefunctions is known 
up to order n. This approach appears particularly promising to compute high-order deriva- 
tives of the total energy with respect to an atomic displacement but, in the formulation of 



Ref. [16], it is of no practical use when the perturbation is an electric field. In fact the formu- 
las contain the change of the eigenvalues of the Hamiltonian due to the perturbation, i.e. a 
quantity which is ill-defined when the perturbation is an electric field and the wavefunctions 
are Bloch states. 

Recently new methods have been introduced in DFT to solve the electronic structure 
problem, mainly to reduce the number of operations necessary for the numerical solution. 
One of these methodsS is based on a Wannier-like representation of the electronic orbitals 
which are constrained to be localized in finite regions of the real space. The localized states 
are in general nonorthonormal and are obtained from a direct minimization of the total 
energy of the system. The method is very convenient to study systems with many atoms 
since the localization of the wavefunctions allows the computation of the total energy with 
a workload proportional to the number of atoms. At the same time, the application of this 
technique to a periodic solid provides a good approximation for the Wannier functions which 
are usually difficult to obtain with other techniques. In Ref. [IB] it was shown that the center 
of these Wannier-like functions yields the correct polarization of the system,0 and that their 
localization property can be conveniently used to study the behavior of a periodic insulating 
solid inside a uniform electric field. This approach allowed the computation of the physical 
properties of a solid under a finite electric field. The derivatives of the energy with respect 
to the electric field were computed by means of accurate finite difference calculations.!!! 



In this paper we further extend the approach of Ref. [18] and Ref. [16| and we give a method 
to compute analytically high-order derivatives of the energy. In particular we give a very 
general derivation of the 2n + 1 theorem, which does not require the definition of a specific 
energy functional, contrary to what was done in previous work. We then apply our result 
to obtain the expressions for the linear and nonlinear susceptibilities in the Wannier-like 



representation of the electronic orbitals.li3tL3 The resulting expression for the second-order 



susceptibility is much simpler than the one obtained by standard perturbation theory!!! 
because the use of the 2n + 1 theorem allows us to express this third-order derivative of the 
energy only as a function of the first order variation of the wavef unctions. Furthermore the 
use in the 2n + 1 theorem of Wannier-like functions instead of Bloch eigenfunctions,^ gives 
an expression for the second order susceptibility which is well defined also for the case when 
the perturbation is a uniform electric field. 

We apply our results to a ID model Hamiltonian to test the convergence properties 
of the proposed algorithm. We compute analytically the first-, second- and third-order 
derivatives of the total energy with respect to a uniform electric field and we compare the 
results with those of the finite difference calculations. The third-order derivative is computed 
for an arbitrary field, so that the fourth-order derivative is available as well through finite 
differences. The simplicity of our method makes it very well suited to compute high-order 
derivatives of the total energy for real materials in the framework of DFT. 

We start with a general proof of the 2n + 1 theorem valid for an arbitrary total energy 
functional ^(w, A), where w is a vector whose elements are the coefficients of all the occu- 
pied wavefunctions on a given basis and A is a parameter measuring the magnitude of the 
perturbation. We restrict ourselves to energy functionals where no explicit constraints — as 
for example those of orthonormalization — are imposed on the components of w.0 For a 
given A the total energy is defined as the minimum of E(w, A) with respect to w. If A is 
varied from A 1 - - 1 to \^ + AA, the vector w which minimizes the energy functional will change 
from w(°) to w^ + Aw. We can expand the total energy around w^ by a Taylor series: 

£ ,w<»» + Aw, A<»> + AA) = gg^ g^M (Aw) > (AA) , (1) 

where we use the notation: |^-(Aw) fc = ( Yli ^ w i~i^) k 'E '■ At a given AA the force is defined 
by: 

_ ggwW + Aw, AW + AA) ~ « 1 ^(wW,AW) _ 
f " p) + Aw) " (fc-l)!p! W (AW) (AA) • (2) 

The vector Aw is the solution of the equation obtained from the extremum condition f = 0. 



We now define f ( n > and as the force and the energy to order n in AA. Explicit expressions 
of these quantities are obtained by writing Aw as: 

Aw = w« + w< 2 > + . . . , (3) 

where is of order (AA) n , and by separating the various orders in Eq. ([!]) and Eq. (0). 
w^ n ) is obtained from the equation f( n ) = 0. Using these definitions the proof of the 2n + 1 
theorem is straightforward. 

Since the term quadratic in is of order (AA) 2Z , £'( 2n+1 ) can contain w^ only at linear 
order if I > n. Under the same condition, we show that the coefficient of w^) in 

E (2n+l) 

is zero. To show this it is useful to single out w^ from the product (Aw)* appearing in 
Eq. ([]]), using the relation: 

(Aw)* = (Aw - w« + w«) fc = k'^Aw)^ 1 + (Aw - w«) fc + o((AA) 2n+1 ), (4) 

which is valid for I > n. The only term which is linear in w^ is the first term of the r.h.s. 
of Eq. (U). Inserting this term in Eq. (p]) and recalling the definition of f, Eq. (0), we can 
write £( 2n+1 ) as: 

E (2n+1) = w (2n+l) f (0) + _ _ _ + w (J)f(2n+l-i) + _ _ _ + w (n+l) f (n) + p(2n+l) ( w (l) ? _ _ _ ^ w (n)^ (5) 

where p( 2n+1 ) is a polynomial of degree 2n + 1. Since fW = for every i, due to the 
extremum condition, we obtain 

E (2n+1) = p(2n+l) (w (l) ; ___ ;W (n) ) _ (g) 

This completes the proof of the 2n + 1 theorem. We note that our formulation does not 
require any hypothesis on the specific form of the energy functional. Therefore our proof 
can be applied to DFT, to correlated wavefunctions, and also in contexts other than quan- 
tum theory. Furthermore the present approach combined with a functional with implicit 
orthonormalization constraints^ can be used to derive the perturbative expansion in cases 
where the standard approach is cumbersome, e.g. the case of DFT when the atoms are 
described by Vanderbilt pseudopotentialsS 
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We now apply the above ideas to the computation of the linear and nonlinear susceptibil- 



ities. As explained in Ref. 18 it is possible to define a total energy functional for a periodic 



insulating solid in a finite electric field as: 

E{\w ),F) = J2( w o\H + eFx\w l ){25 0jl -{w l \w )), (7) 
i 

where H is the unperturbed Hamiltonian of the solid, F is the electric field, x is the position 
operator, e is the electron charge and \wi) are Wannier-like functions associated to the 
direct lattice vector Ri, which are in general non-orthonormal. The Wannier function \wi) 
is obtained by translating the function centered at the origin by a vector Ri. |u> ) is free to 
vary within a real space localization region (LR) of radius R c centered at the origin and it is 
set equal to zero outside LR. For simplicity in Eq. (f7|) we assume that only the lowest band 
is filled, that the system is one-dimensional and that the total energy describes independent 
electrons. Self consistency does not yield any additional problem. We stress here that the 
expectation value of x is well defined for any finite cut-off radius R c . Furthermore we note 
that even if no orthogonality constraints are imposed on the \wi), at the minimum they 
become approximately orthonormal as shown in Ref. [L7|. 

We now recall that the linear and the quadratic susceptibilities x^ an d X^ 2 \ are obtained 
as |xW(A.F) 2 = E^ and |x^ 2 )(AF) 3 = E^ where E^ is the variation of the energy 
functional given in Eq. (|7|) to order n in the perturbing field AF. From Eq. (U) and Eq. (||) 
with AA = AF, we obtain the expressions: 

5*< W 40(w«.>)* + ^w<^, (8, 

where we used the fact that the total energy functional is linear in the electric field. The first 
order variation of the localized orbitals is obtained either from the condition f M = 
as the solution of a linear system:! 
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or from a direct minimization of E^ 2 \ Eq. (§), with respect to as in Ref. [K| All the 
energy derivatives appearing in these formulas are evaluated at the unperturbed vector w^ '. 

We applied the above results to a ID model with Hamiltonian H = — V 2 + V(x) where 
V(x) is a periodic potential with period 3, i.e., V(x + 3) = V(x). We chose V(x) = —A 
if x E (-1.5, -0.5], V{x) = a - A if x E (-0.5, 0.5], and V(x) = if x E (0.5, 1.5]. The 
parameter A is kept fixed at the value A = 4 and a varies between a = and a = A. 
At the two limiting values, the model has inversion symmetry, and therefore x = 0. 
Otherwise the parameter a tunes the value of x ■ We discretized the wavefunctions w{x) 
on a iV-point mesh Xi with equal spacing Ax. In this representation the action of the 
Laplacian operator on the wavefunctions is modelled as a finite difference: V 2 w(xi) = 
(w(xi + i) + w(xi-i) — 2w(xi))/ (Ax) 2 . All the calculations are made with Ax = 1/3. 

In Fig. [l| we show the x values computed from the analytical derivative of the total 
energy, Eq. These are compared with the x values obtained from a numerical differ- 
entiation of the polarization P = jb computed at finite electric fields. The two results are 
in perfect agreement. We note that the finite difference calculation yields the same x^ as 
a perturbative approach based on Bloch wavefunctions. El 

In the same figure we also show a comparison between the values of x^ 2 \ as obtained 
from the analytical derivative, Eq. (||), and from a numerical differentiation of the values of 
X^ computed at finite electric fields. Even in this case the two calculations are in perfect 
agreement. 

In Fig. |^ we present the values of x^ as a function of the electric field for a = 2: its 
linearity around the origin allows us to extract from the slope of the curve the value of 

v (3) _ 2 dx (2) 
a. 3 dF ■ 

All the above calculations have been done with an R r value such that the LR includes 



seven unit cells. As shown in Ref. [18] the main concern in this type of calculations is the 
convergence rate of the studied quantities with respect to the dimensions of the LR. In Fig. |3| 
we show how the values of x^ converge as a function of the size of the LR. 



7 



In conclusion we gave a new proof of the 2n + 1 theorem. We applied it to the total 
energy of an insulator in a uniform electric field where the wavefunctions are described by 
localized orbitals, with no explicit orthonormalization constraints. We provided a method to 
compute analytically the first and second-order susceptibilities, which is much simpler than 
the standard approach. We tested its accuracy and convergence properties in a ID model 
Hamiltonian. We believe that the application of this method to DFT and to an arbitrary 
kind of perturbation could open the way to a simple and reproducible computation of high- 
order derivatives of the total energy. This will be an efficient way to compute important 
properties of real materials such as the Raman tensors or the nonlinear susceptibilities even 
in systems with complex unit cells. 

We gratefully acknowledge many usefull discussions with A. Baldereschi, R. Resta and S. 
Scandolo. We thank R. Car, A. Canning, P. Fernandez, and R. Resta for a critical reading of 
the manuscript. We acknowledge support by the Swiss National Science Foundation under 
Grant No. 21-31144.91. 
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FIGURES 

FIG. 1. Linear (dashed line) and quadratic (solid line) susceptibilities of the model system 
computed analitically with Eqs. (8), (9). The results obtained from numerical differentiation of 
the polarization (solid squares) and of the linear susceptibility (open squares), both computed in 
finite electric fields, are also shown. 

FIG. 2. Quadratic susceptibility as a function of the electric field for a = 2. 

FIG. 3. Quadratic susceptibility as a function of the parameter a for several dimensions of 
the localization regions. The curves refers to a localization region equal to three (long dashed), 
five (dotted), seven (dashed) and nine (solid line) unit cells, respectively. 
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